Epicyclic Oscillations of Fluid Bodies: 
Newtonian Non- Slender Torus 

Omer M. Blaes^ 

blaesOphysics . ucsb . edu 

Eva Sramkova^ 

sram_eva(§centrum . cz 

Marek A. Abramowicz^'^'^ 

marekOf y . Chalmers . se 

Wlodek Kluzniak^'^ 

wlodekScamk . edu . pi 

and 

Ulf Torkelsson^ 
torkelQf y . Chalmers . se 

ABSTRACT 

We study epicyclic oscillations of fluid tori around black holes (in the 
Paczyhski-Wiita potential), and derive exact analytic expressions for their radial 
and vertical eigenfrequencies and to second order accuracy in the width 
of the torus. We prove that pressure effects make the eigenfrequencies smaller 
than those for free particles. However, the particular ratio Vz/^r = 3/2, that 
is important for the theory of high frequency QPOs, occurs when the fluid tori 
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epicyclic frequencies Vr-, are about 15% higher than the ones corresponding to 
free particles. Our results therefore suggest that previous estimates of black hole 
spins from QPOs have produced values that are too high. 

Subject headings: accretion, accretion disks — black hole physics — hydrody- 
namics — X-rays: binaries 



1. Introduction 



The Fourier power density spectra of X-ray va riability in Galact i c black hole X-ray bina- 
ries o ft en reveal pairs of high frequency Q POs (e.g.. lStrohmayerll200ll : lRemillard fc McClintock 
m^. iKluzniak fc Abramowiczl J2001al lb[) suggested that these high frequency QPOs are 
caused by a non-linear resonance between two global modes of oscillations in an accretion 
flow in strong gravity (here we denote these modes by 6r, 6z), and pointed out that the 
observed frequencies are in a commensurate (3:2) ratio. This suggestion was developed by 
them and collaborators into the "QPO resonance rn odel" . The model uses the theory of small 



non-linear oscillations (e.g., iNayfeh fc Mooklll979l ). and attempts to explain many observa- 



tional properties of QPOs in X-ray binaries by deriving them directly from the differential 
equations that describe tw o weakly coupled, no n-linear oscillators (for more information, see 
the collection of reviews in lAbramowica l2005bl ) . 



5f + {ujr)'^6r 
6z + iu},)'^5z 



Xr{5r, 6r, 6z, 6z), 
Xz{6r, 6r, 6z, 6z). 



The resonance model does not address, however, the actual accretion flow structure or the 
specific modes of oscillation, information upon which the detailed form of equations ([T]) 
depend. 

One possibility is that radial pressure gradients set up fluid tori in the accretion flow 
which can support discrete, trapped hydrodynamic modes. That oscillations of such tori 
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RezzoUa et al. 


I2OO3I: see also 


Lee, Abramowicz & Kluzniak 


2004 




Rubio-Herrera & Leel 


20051. and 


Blaes. Arras & Fraeilel 


20061). It is not vet clear 



wheth er such tori provide a realistic model for the accretion flow in the steep power law 
state (IRemillard &: McClintockll2006l ). where high frequency QPOs are observed. Nor is it 
clear whether their hydrodynamic modes of oscillation can exist in the presence of mag- 
netorotational (MRI) turbulence. Nevertheless, pressure supported "inner tori" do appear 
to be an ubiquitous flow feature of nonradiative global simulations of MRI turbulence in 
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accretion flows flHawley &: Balbua 120021 : iDe Villiers. Haw 



such an inner torus is shown in Figure [T] (IMachida et al 



20061 ) 



ev &: Krohkll2003r) . An example of 



If torus-like structures do exist in the steep power law state, global epicyclic oscilla- 
tions of these tori are almost certainly the most robust modes, as their existence is de- 
rived from the properties of the external spacetime, not the internal properties of the torus 
( lAbramowicz et al.ll2006l ). While the existence of these modes is independent of the proper- 
ties of the torus, their actual frequencies and eigenfunctions are not. It is this issue which we 
wish to address in the present paper: how the frequencies of epicyclic modes of thick (non- 
slender) fluid tori differ from the epicyclic frequencies of test particles. As we shall discuss 
later in this paper, this question is of direct relevance for an accurate estimate of the black 
hole spin from the measured QPO frequencies. In order to answer this question, we calculate 
analytically eigenfrequencies and eigenfunctions of the epicyclic modes of nonslender tori up 
to the second order in the torus thickness. 

At first sight, vertical and radial epicyclic modes would seem to be a terrible choice for a 
resonance, as test particle orbits in Kerr spacetime are fully separable and there is therefore 
no nonlinear coupling betwe en these motions. Howev e r, tor i behave like test particles only 
when they are very slender. iKluzniak fc Abramowiczl (120021 ) recognized that for nonslender 
tori the frequencies of the epicyclic modes would be modified by pressure. They derived an 
approximate formula for the epicyclic frequencies, radial Ur and vertical u;^, of fluid tori. 



[UJr 
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and 
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(2) 



Here uJl and cu^ are the radial and vertical epicyclic frequencies for particles, Cso is the 
sound speed at the torus center, and the coefficients A^. and are (not exactly specified) 
function s of the equation of s t ate an d the background gravitational potential. Numerical 
work by iRubio-Herrera fc Led (l2005l ) revealed that > and > 0. In this paper we 
analytically calculate explicit forms of A^ and A^. Another motivation for the present work 
is that it is a necessary step toward deriving an explicit form of equations (1). In the slender 
torus limit, there is no nonlinear coupling of the epicyclic modes, but the pressure corrections 
of equations (2) may give rise to nontrivial couplings. 

For simplicity we model general relativistic e f fects throughout this paper with the 
pseudo- Newtonian potential of iPaczyhski fc Wiital (jl980l ). The mathematics of nonslen- 
der tori is complicated, and a Newtonian calculation is a useful first step before attempting 
the calculation in full Kerr geometry. We will publish an extention of our results to the Kerr 
geometry separately, in O. Straub et al. (2007, in preparation). In any case, the exact an- 
alytic results here will be useful for oscillatory r aode id e ntifica tion in numerically simulated 
Paczyhski-Wiita flows, as already attempted by iBursa I (120061 ) and M. Bursa & M. Machida 
(2007, in preparation). 
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This paper is organized as follows. In Section [2] we briefly review the equilibrium struc- 
ture of tori, citing results that we will need later. In Section [3] we demonstrate that radial 
and vertical epicyclic modes exist for a completely general, baroclinic slender torus. In Sec- 
tion H] we then restrict consideration to polytropic, constant specific angular momentum tori 
and derive the lowest order pressure corrections to the epicyclic mode frequencies (second 
order) and eigenfunctions (first order). We discuss our results and present our conclusions 
in Section [51 



2. Newtonian Slender Torus 

Consider an axially symmetric, inviscid rotating fluid body with toroidal topology in 
equilibrium in an external axially symmetric gravitational field. The flow is stationary and 
its velocity only has an azimuthal component, v = Qrcj). In the paper we use cylindrical 
coordinates {r, 0, z} for all calculations. The gravitational field is described by the potential 
$(r, 2;), which we assume possesses reflection symmetry: $(r, z) = $(r, —z). 

Dynamical equilibrium requires that 

Q^r=^ + V$. (3) 
P 

Here p is the pressure and p is the density. Of particular interest is the circle where the 
pressure has zero gradient. We shall call this the equilibrium circle, as it corresponds to a 
balance between centrifugal and gravitational forces, as one can verify by substituting Vp = 
into equation ([3]). It follows from this equation that the circle lies in the equatorial plane at 
a distance where the rotational velocity Qq and the specific angular momentum £0 of the 
flow have their test particle ( "Keplerian" ) values Qxiro) and ixiro). Let us introduce the 
effective potential of a test particle with specific angular momentum io, U = ^ + il/{2r^). 
The equilibrium circle lies at its minimum and the Euler equation can be rewritten as 

^r = — ^ + VW. (4) 

r-^ p 

The full equilibrium structure can easily be derived from this equation in cases where the 
pressure can be e xpressed as a f unction of density alone (a barotropic or pseudo-barotropic 



equilibrium, e.g. iTassoull Il978l ). In this case, it is possible to find a potential H such 
that VH = Vp/ p. The left-hand side of equation (jl]) can then be expressed as a gradient. 
Moreover, since this gradient has only a radial component, the corresponding potential is a 
function of r only and also the angular momentum ^ must be function of r only - the specific 
angular momentum of the flow is constant on cylinders. 
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We now assume that the equihbrium pressure and density obey a polytropic relation: 



p (X p 



l+l/n 



Let us define two potentials 



H 



dp 
P 



(n + 1 



P 



and \E' = - 



Substituting these into equation (j4]), we obtain the Bernoulli equation in the form 



W + ^ + (ra + 11 



P 



const. 



(5) 



(6) 



The constant can be evaluated by considering the equation at the equilibrium point. Then 
we find 



p 
P 



Po 



nc: 



sO 



Po 



fir,z) 



(7) 



Here c^q = (n + l)po/(^Po)) so that Cso is the adiabatic sound speed Cg evaluated at the 
equilibrium point if the barotropic equilibrium happens to also be isentropic. The pressure 
and density profiles are given by p = po/" and p = po/"'*'^- Surfaces of constant pressure 
and density coincide with surfaces of constant /. 

It is useful to examine the behavior of the function / in the vicinity of the equilibrium 
point. For this purpose let us express the coordinates r and z as r = {1 + x)ro and z = yvQ. 
The equilibrium point corresponds to x = y = 0. For small x and y we have 



/ = 1 



X 



dz"^ 



y 



2io 



de 

dr 



X 



The derivatives of the effective potential come from its expansion. The first derivatives are 
missing because the equilibrium point corresponds to a minimum of the effective potential. 
The mixed second derivatives vanish due to the reflection symmetry. The term containing 
the derivative of the speciflc angular momentum comes from the expansion of the potential 
\E'. The second derivatives of the effective potential with respect to r and z give us radial 
and vertical epicyclic frequencies Ur and Uz in the equilibrium point. Let us express them as 
fractions of the orbital angular velocity Qq at the equilibrium circle, Ur = uir^o, = ^^^^o- 
Then we find 

2ro fdf 
4 \dr/Qj 

where = {2nc^Q) / (rQfll) . The surfaces of constant pressure and density have elliptical or 
hyperbohc cross- sections depending on the sign of the square bracket. In fact, it is possible 
to express the radial epicyclic frequency using the gradient of the test particle ( "Keplerian" ) 
specific angular momentum as ci)^ = (2ro/£o)(c^^K/c^^)o- Therefore, when the specific angular 
momentum of the fiuid increases more slowly than "Keplerian" (with increasing r), the 



/ = 1- 



1 



2,-22 

X +uj,y 



(9) 



- 6 - 



equipressure surfaces have elliptical shape and the equilibrium point corresponds to the 
center of the torus with maximal pressure and density. When the increase is faster the 
equilibrium point is a cusp which corresponds to a saddle point in the pressure and density 
profiles. 

Assume that the specific angular momentum distribution i{r) is such that the equilib- 
rium circle is at the center of the torus. The surface of the torus is the equipressure surface 
where / = 0. Generally it can be far from the center - in the region where our approximation 
(small X and y) is not valid. However, it is clear from equation that our approximation 
is valid in the whole torus if (3 is of the same order as x and y, i.e., when the flow is highly 
supersonic. In this limit the torus becomes infinitesimally slender and its equipressure sur- 
faces have elliptic shape. In the particular case of a constant specific angular momentum 
distribution, these ellipses have semiaxises in the ratio of the epicyclic frequencies. Moreover, 
when the external gravitational field is Newtonian, $ oc l/(r^ + z'^Y^'^, these ellipses become 
circles. It follows that then the torus has a circular cross-section with radius (3rQ. 



3. Epicyclic Modes of Baroclinic Slender Tori 

Keeping the assumption that the torus is made of an ideal fluid, we allow it to have 
arbitrary specific angular momentum i{r, z) and entropy distributions. Surfaces of constant 
pressure and density need not coincide, and rotation need not be constant on cylinders. The 
equilibrium configuration must still satisfy equation (jl]), but now the right hand side of that 
equation cannot necessarily be expressed as a gradient, and we do not bother to attempt to 
solve this equation for a detailed equilibrium structure. However, we still assume that the 
torus equilibrium is slender. 

For simplicity, we restrict consideration to axisymmetric perturbations of this equi- 
librium. The equations describing Eulerian perturbations are those of mass conservation, 
momentum conservation, and adiabatic flow: 

dt ^ p dr dr ' 

^ + iiv.v« = o. (12) 

ot r 

dSv^ _ -1 95p ^ 6pdp 
dt p dz p2 ' 
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and 



d5p 2^^P 



+ 5v • (Vj9 - clvp) = 0. 



(14) 



dt ' dt 



Now, for a slender torus, the derivatives of r in the continuity equation become neghgible, 
and we may write instead 



Hence from now on we may consider all vectors as two dimensional in r and z, which 
may be treated as Cartesian coordinates as far as all vector operations are concerned. We 
consider possible modes in which 6vr and Sv^ are spatially constant. Then after differentiating 
equations ffTTj) and f|T3l) with respect to time, and using equations f|T2|) . f|T^ and f|T5l) to 
ehminate (5w<^, 5p, and 6p, respectively, we obtain 




(15) 




(16) 



For a slender torus, we may write 




) 



(17) 



Hence from the equilibrium condition (jl]), we finally obtain 





The radial and vertical components of this equation give us our modes: 




(19) 



and 




(20) 



4. Behavior of Epicylic Modes for Thick Tori 



We now turn to the behavior of the epicyclic modes for thicker tori. To keep things 
simple, we restrict consideration to polytropic tori with constant specific angular momen- 
tum. This is unlikely to be the most physically relevant case, particularly as existing global 
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simulations of accretion flows with magnetorotational turbulence lead t o near-Keplerian dis- 



tributions of angular momen tum (jPe Villiers. Hawley &: KrolikI l2003l : iMachida et al.l 12006 



Matsumoto fc Machidal 120071 ) . However, the eigenvalue problem is self-adjoint for constant 
specific angular momentum, and the calculations are therefore far more straightforward. 
While there will presumably be quantitative differences for more general angular momentum 
distributions, there will still be pressure corrections to the modes which are likely to be 
comparable to those we calculate here. We discuss this issue further in Section 5. 

Because the equilibrium tori are axisymmetric and stationary, we assume that all per- 
turbations have azimuthal and time dependence oc exp[i{m(j) — out)]. The periodicity of the 
solution in requires integer values of m. 



Papaloizou &: Pringld (119841 ) have shown that the perturbation equations are simplest 
when expressed in terms of the variable 



W 



6p 



p{mVt — uj) 



(21) 



where 5p is the Eulerian perturbation in pressure and VL is the local angular velocity of fluid 
in the equilibrium torus. The Eulerian perturbation in fluid velocity is simply 5w = iVW, 
so that W is directly proportional to the perturbed velocity potential. In terms of W, the 
linear perturbations of the torus are governed by the equation 



ld_ 

r dr 



rf 



dr 



+ 



d_ 

dz 



f 



where 



dz 



r2 ^ ^ I3^rlnl 



p~^w = 0, 



2{n + l)po 



(22) 



(23) 



is our slender torus perturbation parameter, and subscript zero refers to the pressure maxi 
mum of the torus. The equipotential function / is given by 

2 



/ 



where 



U = $(r,;z) + 



2r2' 



(24) 



(25) 



We assume the external gravitational potential $(r, z) is symmetric with respect to reflections 
about the equatorial plane, so all odd 2;— derivatives of U vanish in the equatorial plane. 

The Papaloizou-Pringle equation may be written in abstract operator form as 



LW = -1n{uj - mnfW. 



(26) 
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Here oj = a;/Qo and Q = Q/Qq- Also, 



2~ 



(27) 



is a linear operator which is self-adjoint in the inner product 

<W^\W2>=^J dr j dzrr-'W*W2, (28) 
where the double integral is taken over the torus cross-section defined by / > 0. 

4.1. Slender Torus Limit 

Change variables from r and z to 

Then, in the slender torus limit j3 ^ 0, the equipotential function becomes 

f^'^^l-u^y-u;lf, (30) 

where we use a superscript (0) to denote the slender torus limit. Here LJr and lJz are again the 
radial and vertical epicyclic frequencies at the pressure maximum, scaled with the angular 
velocity Qq- 

1 /d''U\ , o 1 rd^U 



The slender torus limit of the Papaloizou-Pringle equation is 

L(o)|y(o)+2na2W^(°) = 0, (32) 

where ao = co^^^ — m and 



i(c).,.o,^ + „^iL,;.o.|i^„^| (33) 



is the slender torus limit of the hnear operator L. This is still a self-adjoint operator with 
respect to the inner product 

< Wi\W2 >^ J dx J dy{f°Y~^Wi^2, (34) 
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where the integrals are now taken over the domain where /*^°^ > 0. Because of this fact, 
the eigenvalues cTq are all real, and the corresponding eigenfunctions W^^^ form a complete 
orthonormal set in which any regular function defined on the domain f^^^ > may be 
expanded. 



Blaed ( 1l985l ) derived the full set of orthonormal rnodes explicitly for point mass poten 



tials in which Ur 



1. 



Blaes. Arras fc Fragile! (120061 ) have derived the lowest order 



eigenfunctions and eigenfrequencies for the more general case considered here. The simplest 
modes are given in Table 1. We are interested in the behavior of the two epicyclic modes 
( J = 1 and 2 in Table 1) as the torus becomes thicker. In the slender torus limit. 



and W^2°^ = a2ye 



(35) 



These are modes in which the fluid velocity is constant on torus cross-sections and entirely 
radial or vertical in the case of wf^ or w!^, respectively. These modes correspond to radial 
and vertical epicyclic oscillations of a test particle in a circular orbit. The correspondence 
between the fluid modes and test particle epicyclic oscillations can be confirmed by looking 
at the case of a general potential. Substituting equation fl35|) into the Papaloizou-Pringle 
equation fl32l) . we find that they are still solutions of the more general problem. The cor- 
responding eigenfrequencies in the corotating frame {uj — mVLo) are indeed the radial and 
vertical epicyclic frequencies Ur and uj^ in the case of Wr and Wz, respectively. 



4.2. First Order Perturbation Theory 

We now expand the general eigenvalue problem for thick tori in a power series in 13: 

W = + (3W^^^ + (3^W^'^^ + ... (36) 
cu = c^(°)+/3c^«+/3V2) + ... (37) 
/ = /(°)+/5/(^)+/?¥^) + ... (38) 

n = i + /?n« + /52n(2) + ... (39) 

L = L(o)+/3L«+/?2l(2) + ... (40) 

Then if we expand the Papaloizou-Pringle equation (1261) to first order in /?, we obtain 
+ L(°)iy(^) = -2ncxlW^'^ - Anaoito^'^ - mn^'^)W^^\ (41) 
We now expand W^^^ in our orthonormal basis of zeroth order eigenfunctions, 

W('^ = J2b,W^'\ (42) 
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and take the inner product of equation ( HTl) with a particular eigenfunction wf'^ . If J labels 
the particular mode of interest, then this gives us an equation for the first order correction 
to the eigenfrequency: 



UJ 



(1) 



(43) 



If J labels a different mode from the one of interest, then we get an equation for the complex 
coefficients in the first order eigenfunction: 



1 



M(^oj - ^o) 



< 



> 



(44) 



Now, for a constant specific angular momentum torus. 



n 



^0 _rl 



(45) 



Hence, we find Vl^^^ = —2x. The first order correction to the equipotential function is given 
by 

(46) 
(47) 



where 



_ ro fd'U\ - _ro f d^U 



'0 \ / 

Finally, the first order correction to L is given by 

ax^ \ ox J ox oy"^ oy oy 



Hence for the radial epicyclic mode, W^^ = aix, we have 

L(i)|iyf ) >= ai [1 - {oof. + nUrrr) x^ - {^l + nLAr,-^ f] 



(48) 



and 



On the other hand, for the vertical epicyclic mode, w!^"* = a2y, 



and 



(49) 
(50) 

(51) 
(52) 
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In both cases, we immediately see from equation fH3l) that the first order correction to the 
frequency vanishes. 

For the radial epicyclic mode, equations (jll]), fH^ and fISU]) give three nonzero coeffi- 
cients for the first order corrections to the eigenfunction, corresponding to the J = 0, 4, and 
5 eigenmodes of Table 1. Adding these together, the resulting radial epicyclic eigenfunction 
Wi turns out to be given by 



X 



ax 2u}f[i2n + ?>)ul - {n + l)uj\ 

+ \{2u}1 — Uj"^) {—OjI ± SmniDr — nUrrr) ± ^mUJrOjl — ulUrrr + ^lUrzz\ 
+ [ujlUrrr " ^^l^^l T ^mUJr^l - + \)ujllArz^ Ujlf^ + 0{fi'^), (53) 

where a\ is the normalization constant from Table 1. 

The vertical epicyclic mode is much simpler. Equations (jH]), fl5T]) and fl52l) give only 
one nonzero coefficient, corresponding to the J = 3 eigenmode in Table 1. The resulting 
eigenfunction is then given by 

— = y + 4 (±4mu;, - Urzz) xy + 0{[3^). (54) 
a2 ut: 



4.3. Second Order Perturbation Theory 

Because the first order frequency corrections vanish, we have to go to second order to 
determine the finite pressure corrections to the epicyclic mode frequencies. 

The second order expansion terms are OS^^ = Sx"^, 

f^"^^ = — — [x^Urrrr + Qx'^v'^^rrzz + V^^zzzz) ; (55) 



and 



ax^ \ ox J ox oy^^ oy oy 



Here 



"'^^^^^4 v^Jo' ""'^^^4 v^^Jo' '''''^ ^''''^^,\^^), ^^^^ 



The second-order terms in the Papaloizou-Pringle equation fl26|l give 
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-4nao(c^^'^ - mn(2))|^(o) _ 2n{uj^^^ - mn^^^fW^'^\ (58) 

Once again, we expand the second order eigenfunction in terms of the zeroth order 
eigenfunctions, 

W^^) = Y,c,wf\ (59) 

i 

Then on taking the inner product of equation (!58|) with the zeroth order eigenfunction of 
the mode of interest, we find that the second order correction to the frequency is (for the 
case when the first order correction a;^^) vanishes) 

0.(2) = --^(< W^(°)|[L(2) - 4naomn(2) + 2nm\mf]\W'~''^ > 

+ J2bj< -4nmaon(^)]|W^f^ >|. (60) 

j ^ 

For the radial epicychc mode, we find 



2 



8{n + 2)uj!^.lj^ I [[2n + 3)uj^ — [n + 1)UJ:;\ 

X ^^2to^tol[{n + 2)ujI - {2n + 5)cj^] + 2Ujlup-liArrr\{n - \)ujI - (2n - l)cu^] + 4cu^cu^W,,, 

+cu^W,2^J(5 + 6n)cj^ - (3n + l)tu^] + w^W,„.W,^^[(2n - l)w^ - (n + l)cu,2] + 2cu,^(n + \)Ul^^ 
±^mu)lu)l\{ln + 15)cj2 - (I4n + 37)w^] ± 4mcu^cu^W^^^[(7n - \)ujI - (14n + h)u?^ 

±AmiufUrzz[in + 1)0^ - {2n - 5)cu^] + 8m2cu,2cu^[(9 - 7n)LU^ + (14n - ll)cu^] 

+0{(3'). (61) 

The frequency of the vertical epicyclic mode is 



2 



Ul2 = ±i^z + m =F —77 i -Lo'^LolUrrzz - ^t^zzzz - 2u'^U!lUrzz 

8ci;,^ti;f(n + 2) { 

+Urzz[^1i^rrr + {2ojI + 3ci}^)Wrzz] =F 4mc<}2 [cU^CU^ + CU^l^rrr + (4C()^ + 3iu'^)l(rzz] 

+2m'ujl{lQul + 4cu2 - u^)] + 0(/33). (62) 



4.4. Spherically Symmetric Point Mass Potential 

In this case we have 

OJy, — Ll) ^ — 1, lAy^'f — — 6, ^rzz ~ — 
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For the radial epicyclic mode, 



cui = ±1 + m =F 



and 



Wi 



4(n + 2)2 



(53n + 6 ± 7Qnm ^ 8m + 27nm^ - lOm^) + 0{(3^) 



X + 



n + 2 
'3 



—5 =F 4m + (Bra + 1 ± 4m ± 4nm)x 



+ ( -n - 2 =F 4m ) 



For the vertical epicyclic mode, 



a;2 = ±1 + m =F 



and 



4(n + 2) 
2n(ra + l^T^/^ 



(33 ± 52m + 19m^) + 0{P^) 



TT 



[y + (3{3 ± Am)xy] + 0{(3'). 



(63) 



(64) 



(65) 



(66) 



(67) 



It is perhaps interesting to note that for n = 3 (i.e. radiation pressure dominated) tori, 
the axisymmetric (m = 0) radial and vertical epicyclic modes frequencies continue to be 
degenerate to the order of accuracy. 



4.5. Pseudo-Newtonian Potential 



In the case of the iPaczyhski fc Wiital (Il980l ) pseudo- Newtonian potential, in units where 
c = G = M = l, 



we have 



Ur 



- 6 

ro-2' 

ro-2 
4(3rg - 4ro + 2) 

(ro - 2)2 



(r2 + z2)i/2_2' 

6(r2 - Sro + 8) 



(ro - 2)2 ' 
12(3r3 - 30r2 + 60ro - 40) 

3(3ro - 2) 



and 



(68) 



(69) 



These expressions may be substituted into equations flMl) and (1^ to obtain 

rather complicated, explicit expressions for the epicyclic mode eigenfunctions and eigenfre- 
quencies as a function of location of the torus pressure maximum and (3. Rather than present 
those expressions here, we use them to numerically evaluate the mode frequencies. 



- 15 - 



4-5.1. Axisymmetric Modes 

Figure [2] shows the behavior of the axisymmetric {m = 0) radial and vertical epicyclic 
mode frequencies. The /? = curves in these Figures correspond to the test particle frequen- 
cies. For a given pressure maximum radius tq, there is a maximum value of (3 beyond which 
equilibrium tori cannot exist. This limit is indicated by the dashed curves in the Figures. 

Note from Figure [2] that the axisymmetric vertical epicyclic mode frequency for nonslen- 
der tori exhibits a maximum value as a function of ro, in contrast to the behavior of the test 
particle frequency in a pseudo-Newtonian potential. This gives rise to interesting behavior 
of the ratio of the vertical to radial mode frequencies, as shown in Figure [31 Once again, the 
dashed line indicates the limit beyond which equilibrium tori can exist. For small values of 
(3 the ratio of test particle frequencies rises monotonically toward smaller radii, and there is 
therefore a unique radius at which the frequency ratio can take on any specific value. The 
dotted line shows the 3:2 frequency ratio that is found in the high frequency QPOs. As the 
torus thickens, the radius at which the 3:2 commensurability occurs moves inward, and as 
shown in Figure H] (the left part), each of the mode frequencies increases. This continues 
until P = 0.134589, when suddenly tori at two different radii display a 3:2 commensurability 
between the axisymmetric epicyclic modes. The inner torus produces higher frequencies than 
the outer torus. As (3 increases further, these two tori move toward each other, converging 
at /5 = 0.138079. For still thicker tori, there is no radius at which the axisymmetric epicyclic 
modes are in a 3:2 ratio. Hence we conclude that the axisymmetric epicyclic modes can 
represent both observed high frequency QPOs only if the torus is not too thick. The right 
part of Figure H] shows the equipotential surfaces of the torus that exhibits the highest such 
mode frequencies while still retaining the 3:2 commensurability. 

The analytic, poloidal velocity fields of the axisymmetric radial and vertical epicyclic 
modes of this nonslender torus are shown in Figure [51 The radial mode involves some 
vertical expansion on outward displacements (and corresponding vertical compression on 
inward displacements). The reason for this is that the vertical tidal gravity is less at larger 
radii, and so the torus expands under vertical pressure gradients. For the vertical epicyclic 
mode, the motions of at least the outer parts of the torus are also easy to understand. 
The upper half of the torus moves radially outward (and the lower half moves radially 
inward) on upward vertical displacements. This is because the radial component of the 
gravitational field decreases as one moves off the midplane, and so centrifugal and pressure 
gradient accelerations drive outward radial displacements in the upper half of the torus. 
Similarly the increasing gravity as the lower half moves toward the midplane causes inward 
radial displacements. The pressure forces exerted by these radial motions act to oppose the 
vertical motions of the inner parts of the torus. This effect is so strong that the inner parts of 
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the torus shown in the Figure are actually oscillating in the opposite direction to the zeroth 
order vertical epicyclic mode. Our perturbation expansion may therefore not be accurate for 
this thick a torus. It nevertheless makes sense that the vertical displacements in the vertical 
epicyclic mode of a radially extended torus cannot in general all be in the same direction 
simultaneously, as there are no forces that would maintain such radial coherence. 

Increasing the thickness of the torus (increasing (3) always decreases the axisymmetric 
mode frequencies so that they are less than the test particle frequencies. Our complicated 
analytic formulas for the mode frequencies suggest that the physical reasons for the details of 
this decrease must itself be complicated. However, it is not hard to guess physically why the 
mode frequencies should decrease in general. First, unless the radius of the pressure maxi- 
mum is very close to the innermost stable circular orbit, the center of mass of a nonslender 
torus always lies outside the radius of the pressure maximum. At large radii both of the test 
particle frequencies always decrease with radius. Hence the outward shift of the center of 
mass of the torus as it thickens should alone decrease the frequency of the epicyclic modes. 
At small radii, the test particle radial epicyclic frequency increases with radius, and so one 
might expect that the radial epicyclic mode frequency would increase with torus thickness 
(except perhaps when the torus is very close to the innermost stable circular orbit and the 
center of mass then lies inside the pressure maximum). However, in this regime the mode fre- 
quency appears to depend primarily on the nearby presence of the cusp in the equipotential 
surfaces. There is no radial restoring force at the cusp, and as a result the radial epicyclic 
mode frequency of a nonslender torus near the cusp drops dramatically, as shown in Figure 
|5J Finally, even when the pressure maximum of the torus is close to the innermost stable 
circular orbit r < 6.5M so that the center of mass of the torus moves inward on thickening, 
the vertical epiyclic mode is still biased toward the outer parts of the torus because pressure 
forces reduce the vertical oscillations in the inner parts. As a result, the frequency of the 
vertical epicyclic mode still decreases with torus thickness even in this regime. 



4.5.2. Nonaxisymmetric Modes with m = ±1 



Recently, iBursa I (120061 ) suggested a new pair of modes having the 3:2 commensurability, 
that would better satisfy the observational constraints on black hole spins and masses, in 
particular the most recent spin estimates from spectral fi tting to the X-ray spectrum of GRO 



J1655-40 fIShafee et al.ll2006 



Davis. Done, fc Blaesll2006l ). These modes are the axisymmetric 



vertical and the nonaxisymmetric m = —1 radial mode, whose test particle frequencies occur 
in a 3:2 ratio close to the marginally stable orbit. 



The nonaxisymmetric (m = ±1) radial mode frequencies are shown in Figures [H] and [71 
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Again, the curves with (3 = correspond to the frequency of a test particle, and the dashed 
hne denotes the hmit for tori in equihbrium. Figure [6] shows the sum of the axisymmetric 
radial frequency and Qq (i.e., the m = 1 mode) for different torus thicknesses. With increas- 
ing values of f3, the sum of the frequencies slightly decreases. The situation becomes more 
complicated for the m = — 1 mode. The frequency difference between Qq and the axisym- 
metric radial epicyclic frequency first increases with growing /?, but then it decreases with 
increasing (3 in the interval 10.895M < r < 15.573M (Figure [7]). Finally, for radii greater 
than 15.573M, the frequency difference increases with f3 again. 

The ratio of the axisymmetric vertical frequency to the m = — 1 radial mode frequency 
for different torus thicknesses can be seen in the left panel of Figure [H Unlike the resonance 
for the axisymmetric epicyclic modes, the radius at which these two modes are in a 3:2 
ratio moves outward with increasing (3. This corresponds to a decrease of the individual 
frequencies, as can be seen in the right part of the Figure, and the frequencies eventually 
become zero for P = Pmax = 0.725108. For tori with P < Pmax, there exists exactly one torus 
pressure maximum location that generates a 3:2 frequency ratio, and no such location exists 
for still thicker tori corresponding to (3 > Pmax- This would suggest that even fairly large 
tori can exhibit a 3:2 resonance, but the perturbative method is valid only for small f3, and 
any firm results on very extended tori must be obtained by other methods. 



4.6. Consequences for Black Hole Spin Estimates 



The first black hole spin estimates from fitting the axisymmet ric epicyclic mode frequen- 



cies to the obser ved QPQs from GR Q J1655-40 were reported by lAbramowicz fc Kluzniak 



( I2OOII ). Recently iTorok et al.l (120051 ) have used the orbital resonance model to estimate the 
spins of all three microquasars with known masses. The most tightly constrained spin was 
in GRO J1655-40, where a/M was found to lie between 0.93 and 0.99. This spin value that 
was derived from a resonance between the axisymmetric epicyclic modes is not in accord 
with the independent spin measurement s, determined by fitting the X-ray spectral continua: 
a/M ~ 0.65 — 0.75 (IShafee et al.ll2006l ). a/M ~ 0.62 or less if the disk inclination is free 



( Davis. Done, fc Blaea 



20061 ) for GRO J1655-40. 



However, the resonance model estimates were based on the epicyclic frequencies for free 
test particles. Figure H] shows that the axisymmetric epicyclic frequencies at the resonant 
radius will be higher for a nonslender torus than for a test particle. The maximal frequency 
shift due to pressure effects is around 15 percent of the test particle frequency. The previous 
studies have therefore most probably overestimated the actual value of the spin of the black 
hole, which explains some of the discrepancies with other methods. 
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On the other hand a resonance between the axisymmetric vertical epicychc mo de and 



thenonaxisymmetric m = —1 radial mode yields a black hole spin for GRO J1655-40 ((Bursa 



20061 ) compatible with other recent spin estimates although Bursa used the frequencies for 
free particles. As seen from Figure El the frequencies of these two modes at the resonant 
radius decrease with growing torus thickness, contrary to the axisymmetric modes. Therefore 
in this case the black hole spins should be higher than previously estimated. The maximal 
shift of the frequencies (and consequently the black hole spin) can be very large, but our 
perturbative method is valid for small f3, and is not reliable for (3 ^ 0.7. 

The above modifications of the spin estimates are based on purely Newtonian calcula- 
tions using the Paczyhski-Wiita gravitational potential. In order to explore the real behavior 
of thick tori orbiting around rotating black holes, one needs to solve the hydrodynamical 
equations in the Kerr metric. 



4.7. Comparison with Numerical Simulation 

We compare our analytically calculated frequencies for the axisymmetric epicyclic modes 
with those from hydrodynamic numeric al simulations. The nu merical simulations were 



performed using the ZEUS-2D code by IStone &: NormanI (119921 ). The radial, and verti- 



cal epicyclic oscillations, respectively, were excited by applying a purely radial, or vertical 
constant velocity field perturbation, respectively, to the torus at t=0. Then we Fourier an- 
alyzed the radial and vertical positions of the center of mass in order to find the oscillatory 
frequencies. 

Figure M compares the results of our analytic expressions (1641) and (I66p with m = 
to the epicyclic mode frequencies measured from the simulations as a function of torus 
thickness. The agreement is excellent for (3 ^ 0.2. Figure [TO] similarly compares the analytic 
m = frequencies from our perturbation theory to the simulation results for tori in pseudo- 
Newtonian potentials. Once again, there is good agreement for f3 < .2. The general trend 



that t he frequencies decrease with increasing f3 was also observed by iRubio-Herrera fc Lee 



(120051 ) in their simulations. 



5. Discussion and Conclusions 

In this paper we have assumed that (possibly nonaxisymmetric) vertical and radial 
epicyclic modes of fluid tori are responsible for the commensurate pairs of high frequency 
QPOs in black hole X-ray binaries. We have derived exact analytic formulae for the mode 
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eigenfunctions and frequencies. We find that the frequencies are always below those corre- 
sponding to free test particles, an effect that must be taken into account before identifying 
such modes with observed QPOs. 

While our analytic calculations are exact, they apply to very simphfied and idealized 
configurations: 

(1) Our calculations were restricted to Newtonian mechanics. This is a necessary first 
step to doing the fully general relativistic calculation in Kerr spacetime (Straub et al. 2007, 
in preparation). Moreover, Newtonian potentials that mock u p various aspects of the physics 



of Schwarzschild spacetime are used by many numericists (e.g. iLee. Abramowicz fc Kluzniak 



2004 : iMachida et al.ll2006l ). and our results may be applied directly to their simulations. 



(2) We considered isolated, non-accreting, non-magnetized, polytropic tori. The flow 
configuration of the accretion disks in X-ray binaries when they exhibit QPOs is still far from 
clear, but pressure-supported torus-like configurations remain an interesting possibility. Such 



confi g urations are seen as " i nner tori" in global MRI simul ations (jPe Villiers. Hawley fc Krolik 



20031 : iMachida et al.l l2006l : iMatsumoto fc Machidal 120071 ) . These inner tori have magnetic 



energy densities much less than the gas internal energy density, so magnetic fields are less 
important to the overall hydrostatic structure than gas pressure. (Strongly magnetized con- 
figurations may however exist in real systems, and our treatment in this paper would not 
be valid for them.) Among all the possible modes of an isolated torus, the epicyclic modes 
will be most robust to boundary conditions, as they correspond to physical displacements of 
the entire torus. We therefore believe that our calculations for the corrections to the mode 
frequencies will be reasonably robust to these boundary conditions. The entropy distribution 
of these tori is completely unknown, and existing simulations tell us little as they do not 
have consistent thermodynamics. Nevertheless, the luminous accretion flows exhibiting high 
frequency QPOs are likely to be radiation pressure supported, and a polytropic configuration 
with n = 3 may therefore be fairly close to reality. 

(3) For mathematical reasons, we also restricted ourselves to constant specific angular 
momentum tori. Such angular momentum distributions appear unlikely to exist in nature, as 
global simulations of disks with MRI turbulence produce inner tori with angular momentum 
distributions which are closer to the test particle ( "Keplerian" ) limit. (Note, however, that 
these simulations currently lack optically thick radiative cooling which is likely to determine 
the overall pressure support of the flow, and therefore departures from Keplerian angular mo- 
mentum.) The Papaloizou-Pringle equation with non-constant specific angular momentum 
is not a self-adjoint eigenvalue problem, and so we do not have a complete set of orthogo- 
nal basis eigenfunctions in order to expand the perturbed modes around the slender torus 
limit. The perturbation equations using the Lagrangian displacement vector as dependent 
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variable can be used to generate a complet e set of orthogonal eigenf unctions, at least within 
Newtonian mechanics (jSchenk et al.ll2005l ). and this may enable an analytic calculation of 
the corrections to the mode frequencies for non-constant specific angular momentum tori. 
We note that fully general relativistic simulations in Kerr spacetime of non-constant specific 
angular momentum oscillating tori exhibit axisymmetric radial and vertical epicyclic mode 
frequencies which are below the test particle frequencies (Fragile 2006, private communica- 
tion), in agreement with the qualitative results we have found here. 

We view the results of the calculations presented here as a necessary tool for identifying 
the modes in fully 3D, time dependent numerical simulations o f magnet ohyd ro dynamic tori. 
A preliminary example of what can be achieved is provided by iBursa I (120061 ) . who used our 
results to iden tify the (m = 1) radi al and vertical epicyclic modes in the recent numerical 



simulations by lMachida et al.l (120061 ). Once identified, our analytic expressions for the mode 
eigenfunctions may also be useful in trying to explain why such modes are excited in simu- 
lations. They may also be u seful in making pre dictio ns for the observed X-ray mod ulation. 
As explored numerically by iBursa et al.l (120041 ) and ISchnittman fc RezzoUal (120061 ) , gravi- 
tational redshifts, Doppler shifts, and lensing of X-ray photons emitted by tori oscillating 
in both the vertical and radial epicyclic modes can produce detectable oscillations in the 
predicted X-ray fluxes measured by an observer. 
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Fi g. 1. — A press u re sup port ed inner torus is appare n t in t he MHD simulations described 
by iMachida et al.l (120061 ) and iMatsumoto fc Machidal (120071 ) . The relative density is color 
coded. The superimposed elliptical shape (a white broken line) corresponds to the analytic 
nonslender torus studied in this paper. 
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Fig. 2. — Left: Frequency of the axisymmetric (m = 0) radial epicyclic mode, i^i = a;i/(27r), 
for an n = 3 torus in a pseudo-Newtonian potential with black hole mass M — IOMq, as 
a function of radius of the torus pressure maximum in units of the gravitational radius. 
Different solid curves correspond to different values of the pressure parameter /3, and are 
labeled as such. Models below the dashed curve are unphysical, as the corresponding tori 
extend outside the critical cusp equipotential. Right: Same as the left Figure, but for the 
axisymmetric vertical epicyclic mode. Once again, physical models lie above (to the right) 
of the dashed line. 
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Fig. 3. — Ratio ^2/^1 = 1^2/^1 of the vertical to radial axisymmetric epicyclic mode frequen- 
cies for an n = 3 torus in a pseudo-Newtonian potential, for different values of f3 and as a 
function of radius of the torus pressure maximum. Again, physical models lie above and to 
the right of the dashed hne. The dotted hne indicates a 3:2 commensurability between the 
two modes. For a slender torus with a particular value of /3 < 0.134589, there is one radius 
for which the axisymmetric mode frequency ratio will be 3:2. For somewhat thicker tori with 
values of /3 satisfying 0.134589 < (3 < 0.138079 = /3max, there are two locations of the torus 
pressure maximum where the mode frequencies are in a 3:2 ratio. No such radii exist for still 
thicker tori with j3 > Pmax- 
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Fig. 4. — Left Axisymmetric vertical epicyclic mode frequency V2 = 1^2/(277) for tori which 
have the axisymmetric vertical and radial modes in a 3:2 ratio, plotted as a function of f3. 
Again, the tori are assumed to have n = 3 and orbit in a pseudo-Newtonian potential with 
M = IOMq. The /3 = 0.134589 point in the Figure corresponds to the intersection of the 
dashed and dotted lines in Figure [31 Right: Meridional cross- sect ion of a nonslender torus 
with j3 = 0.134589 and a cusp at r = 5.0485M. The torus is centered at r = 7.293M where 
the axisymmetric epicyclic mode frequencies ui and 002 for this (3 are in a 3:2 ratio. This 
case corresponds to the crossing point of the dashed and dotted lines in Figure |3l 



Fig. 5. — Poloidal velocity field for the radial epicyclic mode (left) and the vertical epicyclic 
mode (right) of a nonslender n = 3 torus with (3 = 0.134589 and pressure maximum at 
r = 7.293M (the same torus illustrated in the right hand panel of Fig. H]). 
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Fig. 6. — The sum of the axisymmetric radial epicychc mode frequency ooi and the orbital 
frequency Qq at the pressure maximum, for an n = 3 pseudo- Newtonian torus around a 
M = IOMq black hole, as a function of tq for different values of (3. Physical tori exist only 
to the right of and above the dashed line. 
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Fig. 7. — The difference between the orbital frequency VLq at the pressure maximum tq and 
the axisymmetric radial epicyclic mode frequency for an n = 3 pseudo-Newtonian torus 
around a M = IOM0 black hole, as a function of tq for different values of (3. Between the 
radii r = 6M and r = 10.895M the frequency difference increases with increasing /3, and 
physical tori exist only below the dashed line. From r = 10.985M to r = 15.573M, the 
frequency difference decreases slightly with increasing /3, and physical tori exist only above 
the dashed line. For radii greater than 15.573M the behavior is similar to that of the first 
interval. 
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Fig. 8. — Left: The ratio cj2/(^o ^ ^^i) of the axisymmetric vertical mode frequency to the 
frequency of the (m = — 1) radial mode for an n = 3 torus in a pseudo-Newtonian potential, 
plotted for different values of /3 as a function of radius of the torus pressure maximum. 
Physical models he above the dashed hne. For all thick tori having /3 < P^^ax — 0.725108, 
there is one possible location of the torus pressure maximum at which the modes are in 
a 3:2 ratio. For tori with (3 > Pmax, no such pressure maximum radii can exist. Right: 
The axisymmetric vertical epicyclic mode frequency z/2 = cj2/(27r) (for the same tori as in 
the left panel) at the radius where 002/(^^0 — uji) = 3/2, plotted as a function of At 
P = Pmax = 0.725108, the frequency 1/2 goes to zero. 
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Fig. 9. — Left: Comparison between the analytically (1) and numerically (2) calculated 
axisymmetric radial epicyclic mode frequency for an n = 3 torus orbiting in a spherically 
symmetric point mass potential, plotted as a function of (3. Right: Same as the left panel, 
but for the axisymmetric vertical epicyclic mode. 



-30- 




0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 



Fig. 10. — Left Comparison between the analytically (1) and numerically (2) calculated 
axisymmetric radial epicyclic mode frequency for an n = 3 torus orbiting in a pseudo- 
Newtonian potential, plotted as a function of /?. The pressure maximum of the torus is at 
9.2M. For this radius, the maximum value of (3 for which an equilibrium torus can exist, is 
0.492968. Right: Same as the left panel, but for the axisymmetric vertical epicylic mode. 



Table 1. Simplest Modes of the General, Constant Specific Angular Momentum Slender Torus. 



Eigenfunction'^ 



ao 

X bj^ aix 

2 qI a2y 

3 i^r+'^z °-z^y 

4 (2n+l)(.I,;+a.;)-[4n(n + l)(^^-^;)V(a.; + ^;)^]l/^ T "2o_^-2^2 _ ,-2,-2^ _ 2(n+l).D^ 



is an arbitrary nonnegative integer in dex which we are using to lab el the modes. A more physically descriptive 
set of labels can be found in the paper bv lBlaes. Arras &: FragiTg 1 I2OO6I ). 

''The constants {ao, ai, 05} are given in Table 2, and are chosen such that the eigenfunctions are normalized 
in the inner product l|34|l . 
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Table 2. Normalization Constants for the Eigenmodes of Table 1. 
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(1155^)1/2 

1 ao[2(n+l)tD2]l/2 

2 ao[2(n+ l)a;2]l/2 

3 ao[4(n+ l)(n + 2)(I.2<l>2]i/2 
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